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Abstract 

We study the possibility to constrain deviations from Lorentz invariance in dark 
matter (DM) with cosmological observations. Breaking of Lorentz invariance generi- 
cahy introduces new hght gravitational degrees of freedom, which we represent through 
a dynamical timelike vector field. If DM does not obey Lorentz invariance, it couples 
to this vector field. We find that this coupling affects the inertial mass of small DM 
halos which no longer satisfy the equivalence principle. For large enough lumps of 
DM we identify a (chameleon) mechanism that restores the inertial mass to its stan- 
dard value. As a consequence, the dynamics of gravitational clustering are modified. 
Two prominent effects are a scale dependent enhancement in the growth of large scale 
structure and a scale dependent bias between DM and baryon density perturbations. 
The comparison with the measured linear matter power spectrum in principle allows 
to bound the departure from Lorentz invariance of DM at the per cent level. 



1 Introduction 

Elucidating the nature of dark matter (DM) stands as a major challenge of modern cosmology 
and particle physics. One of the basic properties always assumed is that DM satisfies Lorentz 
invariance (LI). In this work we want to analyze how this assumption can be verified from 
the study of cosmological perturbations. 

Lorentz invariance is one of the best tested symmetries of the Standard Model of particle 
physics [IJ. Thus, it is tempting to postulate that it is a fundamental property of all fields 
of Nature including gravity and the dark sectors of the Universe, i.e DM and dark energy. 
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This assumption is very powerful, but also restrictive. When applied to gravitation it has 
far-reaching conceptual implications, like the uniqueness of General Relativity (GR) as the 
LI theory of gravity, see e.g. [21 El IH IS]- However, beautiful as it is, GR suffers from the 
problem of non-renormalizability precluding its interpretation as a UV complete quantum 
theory. This is the essence of the notorious problem of quantum gravity. 

It is conceivable that the eventual theory of quantum gravity will involve violation of 
LI (Lorentz violation, LV, for short) in some form. For instance, it was recently suggested 
by P. Hofava that a UV completion of GR may be possible within perturbative quantum 
field theory at the cost of abandoning LI at very high energies [6]. This approach appears 
promising as it allows for an immediate improvement compared to GR coming from the 
softening of the gravitational amplitudes in the UV0 m. If deviations from LI are present in 
quantum gravity, it is mandatory to understand which are the consequences for the rest of 
fields in Nature 

A first relevant observation is that even when introduced at very high energies, LV has 
also consequences at low energies [10] . Indeed, LV corresponds to the existence of a preferred 
reference frame at every point of space-time. Additional local fields are required to set up 
this frame; those have typically massless excitations and affect physics at any scale. The 
situation at low-energies is encapsulated by the so-called Einstein-aether model [TU |T2] where 
the preferred frame is determined by a vector field with unit norm (aether). This model is 
an effective theory with a cutoff that can be as high as the Planck mass [13] . Restricting the 
aether to its longitudinal component one obtains the khronometric model that was introduced 
in pL4j as the low-energy limit of the Hofava's proposal we just mentioned [TJ [15] (more 
precisely, of its healthy extended version [15]). In both Einstein-aether and khronometric 
cases, the interaction of the aether with the metric field has consequences for gravitation at 
all distances, which leads to observable astrophysical and cosmological effects [121 HZl HE] • 

Precision tests of LI in the Standard Model of particle physics put extremely tight con- 
straints on the coupling of the aether to ordinary matter [H [8]. If similar bounds held for 
the couplings of the aether to other sectors of the Universe, LV would have no effect on 
cosmology. In this paper we adopt the viewpoint that there is a mechanism that enforces LI 
of the Standard Model with the required precision while allowing for sizable LV in gravita- 
tional physics, DM or dark energy. Finding such mechanism while avoiding fine-tuning is an 
important challenge for the model and a few possibilities have been proposed in [191 l20l [21]. 
Some consequences of LV for dark energy were unveiled in [22] , where it was shown that the 
presence of the aether vector field allows to attribute the current acceleration of the Universe 
to a renormalizable operator not sensitive to UV corrections. The present work is devoted 
to a systematic investigation of LV effects in the cosmological evolution of DM. To isolate 
these effects we will assume throughout the paper that the dark energy is represented by the 
cosmological constant. 

We consider a general set-up, with DM described in terms of an effective fluid that in- 
teracts with the aether to account for LV. The simplest DM candidate described within 

^One should keep in mind though that the renormahzabihty of gravity in the strict sense along these lines 
has not yet been demonstrated due to the complexity of the resulting theory. 
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our formalism is the theory of massive non-relativistic particles, but other possibilities, such 
as axionic DM [23], are also covered. We will assume that the velocities of DM are non- 
relativistic during the observable evolution of the Universe. One might think that this would 
suppress all LV effects. However, this is not the case. We will see that the coupling to the 
aether modifies the inertial mass of the DM particles, but does not affect their gravitational 
mass. As a consequence, DM violates the equivalence principle and the dynamics of cos- 
mological perturbations change. The situation resembles that in LI theories with additional 
long-range interaction in the DM sector fM\ 123] . However, the specific signatures of LV allow 
to distinguish between these two cases. 

It is worth pointing out that a future direct detection of DM may yield strong constraints 
on LV in this sector. Indeed, a direct detection would imply a relatively strong coupling 
between the DM and visible sectors. In that case LV would be transferred from DM to 
ordinary matter via loop diagrams and thus would be subject to the tight existing bounds 
on LV in the Standard Model. Similar arguments were used in [261 EZ] to constrain violation 
of the equivalence principle in the DM sector in the LI context. However, in the absence of 
any direct detection so far and since the loop calculations are rather model-dependent, we find 
it useful to study the bounds on LV in DM following from its gravitational manifestations, 
in particular in cosmology. 

Our work is organized as follows. In Sec. [2] we briefly review the Einstein-aether and 
khronometric models. To understand the qualitative effects of LV in DM we study in Sec. [3] 
the dynamics of massive point particles coupled to the aether, concentrating on the Newto- 
nian limit and Jeans instability. We next proceed to the systematic treatment of LV DM 
using the effective fluid description. Sec. H] contains the description of the setup, analysis of 
the homogeneous cosmology and the equations for the linearized perturbations. Sec. con- 
tains the analytic study of the perturbations in various regimes and the qualitative discussion 
of the observational signatures. The results of the numerical integration of the linearized 
equations in a simplified cosmological model are presented in Sec. |6l Section [7] contains the 
summary and discussion of our results. Some technical details of our numerical procedure 
are described in the Appendix. 

2 Gravity with a preferred frame 

To describe LV we assume that at every point of space-time there is a time-like vector Uj, 
that following Ref . ^Jj we will call "aether" . The vector is constrained to have unit norno. 




Thus it does not vanish anywhere and sets the preferred time-direction at every point of 
space-time. This breaks the local Lorentz group SO {3, 1) of GR to the local subgroup SO (3) 

^We use the metric with (+, — , — , — ) signature. Latin indices from the middle of the alphabet take the 
values i,j = 1,2,3, while Greek letters denote the space-time indices. The latter are manipulated with 5^^. 
Objects in bold face are three- vectors. We use units where the speed of propagation of light is c = 1. 




(1) 
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of purely spatial rotations that leave invariant. This way, the introduction of allows us 
to describe LV effects with an action invariant under arbitrary coordinate transformations. 
We will be interested in the dynamics of at large distances, which according to the rules of 
effective field theory is governed by operators with the smallest number of derivatives. It is 
straightforward to see that it is impossible to construct any contribution to the Lagrangian 
with one or no derivatives. Thus LV at low-energies is governed by the action of the Einstein- 
aether model pT| [12], the most general action containing up to two derivatives of m^. 



i? + ir^^^vxv.w" + /(wX-i) > (2) 



where 

ir'^^, ^ cg'^'^g^, + C25i^5; + c,5>;5l + cu^'n'^g^p, (3) 

and the last term with the Lagrange multiplier I has been added to enforce the unit-norm 
constraint. We have included in the above action the Einstein-Hilbert term for the metric 
g^y. The parameter Mq is related to the Planck mass, cf. ([2]), while the dimensionless 
constants Ca, a = 1,2,3,4, characterize the strength of the interaction of the aether with 
gravity. As discussed below, observations require these constants to be much less than 1, so 
we will assume 

\Ca\ < 1 (4) 

throughout the paper. Let us stress again that while we are interested in describing LV, 
the action ([2]) is explicitly generally covariant. This stems from the assumption that the 
non-vanishing aether field represents the only source of LV. In a certain sense, the picture 
is analogous to the spontaneous symmetry breaking of internal symmetric^. Finally, as an 
effective theory Einstein- aether has a cutoff of order Mo|ca|^/^ [13]. If are not extremely 
small this cutoff is only a few orders of magnitude below the Planck mass. 

A variant of the Einstein-aether model is obtained by restricting the aether to be hyper- 
surface-orthogonal , 

^yg^"'^^a^,ya 

In this case the unit-norm constraint is identically satisfied and the Lagrange multiplier term 
in the action ([2]) can be omitted. The expression ([5]) is invariant under reparameterizations 

a I— 7- cr(o-) , (6) 

where a (a) is an arbitrary monotonic function. The scalar field a is assumed to have a time- 
like gradient at every point of space-time. This defines a preferred time-coordinate, hence 



•^There is an important difference though. The condition ([Ij excludes from consideration any LI states 
thus precluding the symmetry restoration. The attempt to make the length of dynamical to allow for the 
restored symmetry phase leads to appearance of a ghost degree of freedom that spoils the theory [28j . The 
absence of a healthy LI phase is common to gravitational models with LV, which raises doubts about the 
possibility to UV complete them within a LI setting, see the discussion of this issue in [521 130] • 
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a is called "khronon" (from the Greek word for "time"), and the class of models including 
the metric and a "khronometric" models. When restricted to the form ([S]), the curl of the 
aether vanishes 

= e/^-V^^VAUp = . (7) 

As a consequence, the four combinations with derivatives of in the action ([2]) are not inde- 
pendent. It is common to eliminate the Ci-term in favor of the other three. The coefficients 
for the last three operators in ^ then become 

A = C2 , /3 = C3 + ci , a = C4 + ci . (8) 

This model naturally arises as the low-energy limit of Hofava gravity [151 [H]. other words, 
Hofava gravity can provide a UV completion for the khronometric theories which potentially 
improves their UV behavior as compared to GR. 

The main difference between Einstein-aether and khronometric models is the number of 
degrees of freedom. Einstein-aether describes three types of massless propagating modes: 
the standard transverse-traceless tensor modes of the metric and the vector and scalar po- 
larizations coming from the transverse and longitudinal fluctuations of the aether. For the 
khronometric case, the transverse vector polarization is absent. In the scalar and tensor 
sectors the two models are almost equivalent0. As in this paper we are mostly interested in 
the scalar cosmological perturbations, we will often use without loss of generality the termi- 
nology of the khronometric model and, in particular, the constants ([H]). One can show that 
the parameters Ca (or ([8]) for the khronometric case) can be chosen such that all modes are 
stable and have positive energy [121 US]- To avoid gravitational and aether Cherenkov losses 
by high-energy cosmic rays [28], one should also require that the velocities of the graviton 
and aether modes are not less than the speed of high-energy particles comprising the cosmic 
rays, which is close to one. 

The phenomenology of theories with action ([2|) has been extensively studied [121 [HI HH] • 
The high precision with which LI is tested within the Standard Model excludes the direct 
interaction of the aether with baryonic matteiH, meaning that the latter couples universally 
to the metric g^j^y. Still, the aether affects the gravitational field of matter sources. At 
the Newtonian level the modifications amount to an unobservable renormalization of the 
gravitational constant that now reads [321 |T6] , 

The deviations from GR for Solar System physics are encoded in two post-Newtonian pa- 
rameters af^^ and af'^^. These have been calculated for the generic aether in Ref. [33] 
and for the khronometric model in Refs. [HI [18]. Observations yield the constraints [31] 

l«r"^l<10-\ |a^"'|<4■10-^ (10) 

The only difference is the presence of an instantaneous mode in the khronometric model, which is absent 
in Einstein-aether 31 . This difference is not relevant for the local physics that we study in the present work. 
^Examples of mechanisms to suppress such interaction are discussed in [TOl [^ . 
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Those constraints are trivially satisfied in GR, where af = 02 — 0- For LV theories 
with generic parameters they imply the condition 

\Ca\ < 10-' . (11) 

This bound is relaxed for certain relations between the parameters. In the generic aether 
model one can impose two restrictions on to make both af and vanish [12] . Then 

one is left with a two-parameter family of theories that are indistinguishable from GR at the 
post-Newtonian level (i.e. for Solar System tests [35]). Remarkably, for the khronometric 
case the same is achieved by the single condition a = 2/3, which leaves the parameters /3 and 
A arbitrary. Further bounds of order 

\Ca\ < 10-^ (12) 

follow from considerations of Big Bang Nucleosynthesis (BBN) [32] and emission of gravita- 
tional waves by binary systems [3S1 [IS] ■ 

To sum up, the theory ([2]) with universal coupling of ordinary matter to the metric g^^ 
is phenomenologically viable provided the stability constraints are satisfied and (jlj) holds 
at the level of (fTT]) or (fT2!) (in the latter case the relations ensuring vanishing of the PPN 
parameters must be fulfilled). 



3 Lorentz violating dark matter: point-particles 

To grasp the physical consequences of LV in the DM sector we first study how the coupling 
to the aether affects the dynamics of gravitationally interacting point particles. This will 
allow us to understand the Newtonian limit of the theory and the effects of LV on the Jeans 
instability. 

In the presence of the aether the relativistic action for a massive point particle can be 
generalized to 

Spp = -mjds F{uX) , (13) 

where 



ds = ^/g^udx/'dx" (14) 
is the proper length along the trajectory of the particle and 

is the particle's four- velocity. F is an arbitrary positive function that we normalize to 
F{1) = 1; the GR limit corresponds to F = 1. Note that particles described by the action 
(|T3|) violate the equivalence principle and, actually, do not follow geodesies of any metric. 
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It is instructive to work out the relation between the covariant description based on the 
action fll3p and the approach commonly adopted in the study of LV in microphysics. In the 
latter approach LV manifests itself for free particles in the modified dispersion relation 



= £\p') , (16) 

where E and p are the energy and momentum of the particle in the preferred frame where 
the aether is aligned with the time direction, uq = 1, Ui = 0, and is a function that takes 
the form, 

c2 2,2 

if LI is preserved. The preferred frame singles out a preferred time coordinate t and through 
( 1T5]) defines the Lagrangian 

Spp = dt Lpp , (17) 



where 

^-"-'"^^Ktt^)' '''' 

The energy and momentum are now related to the velocity by the standard formulas 

For a given function F these equations yield an implicit relation between E and Pi, i.e. the 
dispersion relation of the form f|T6l) . Vice versa, starting from the dispersion relation f|T6|) 
one obtains the equation for the Lagrangian 

i dLpp _ 2fdLppdLpp\ 

We conclude that there is a one-to-one correspondence between the functions F of the co- 
variant approach and the functions S"^ appearing in the free-particle dispersion relation. 
Let us solve Eq. (1201) in the simplest case of the quadratic dispersion, 

^' = m2 + (l + 0p', (21) 

where ^ is a dimensionless constant. Under fairly general assumptions this type of dispersion 
relations with ^ ^ describes the leading effect of LV at relatively low energies |3B] • In this 
case we observe that by the redefinition = V^{1 + equation (120|) reduces to the 

standard relativistic form. Thus we obtain, 




Lpp = -m V 1 - V2 = -m^/ 1 - ^ , (22) 
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and finally 

S„=--/d.^/i±||f?. (23) 

As discussed above, this case encompasses a wide class of physically interesting situations. 
However, since imposing the specific form fl2^ does not simplify the analysis, and also to be 
completely general, we will continue to work with an arbitrary function F. 



3.1 Newtonian limit 

We will assume that DM is non-relativistic during the relevant stages of cosmological evo- 
lution, i.e. it moves slowly with respect to the cosmic frame. The latter is defined as the 
frame where the cosmic microwave background (CMB) is approximately isotropic. We will 
assume that it coincides with the preferred frame set by the background value of the aetheiEI. 
Thus one can use the Newtonian limit to describe the dynamics at subhorizon scales. This 
corresponds to expanding the action f lT^ to quadratic order in the particle three-velocities 
quadratic order in the spatial component of the aethei0 and to linear order in the 
Newton potential 0. The latter appears in the standard Newtonian limit of the metrical, 

^00 = 1 + 20 , ^oi = , g^, = -6ij{l - 2^) , (24) 

where for the moment we have neglected the cosmological expansion. We obtain 



Spp = m I dt 



Y- 



(25) 



where we have denoted 

Y = F'{1) (26) 

and omitted the constant term corresponding to the rest-mass. To understand the effect of 
LV, we consider first the case when the aether fluctuations are negligible, = 0. In this 
case the last term in fl25l) renormalizes the particle's inertial mass 

m I— )■ m(l — Y) . 

On the other hand, the gravitational mass (the source of in fl25|) ) remains equal to m, 
which clearly violates the equivalence principle. To guarantee the positivity of the kinetic 
energy we impose the restrictions m > and y < 1. 



^This is justified since it has been shown [3S] that in a homogeneous expanding universe the aether field 
tends to align with the time direction. 

^ We assume that is of the same order as , as will be verified by the calculation for the self-gravitating 
situation given below. The presence of the direct coupling between DM particles and the aether is crucial 
for that. In the absence of such coupling, as it happens for the field produced by visible matter [331IIH], the 
aether has the order 0{V^) and the dynamics of the system is different. 

^ Other metric perturbations are of higher post-Newtonian order. 
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Let us now consider the generic situation in the Newtonian hmit for a dense medium 
composed of DM particles. It is convenient to introduce the mass density, 

p(x,t) = m^5(3)(x-x^(t)), (27) 

A 

where Xyi(t) is the trajectory of the A-th particle, and rewrite the action (125|) in the form. 



Spp = d^x p 



Y- 



(28) 

One observes that inside the medium the aether perturbations acquire a quadratic potential 
with the central value set by the velocity of the medium. Not to destabilize the aether 
the potential must be positiv^, Y > 0. We can anticipate that due to this potential the 
aether tends to align with the velocity of the medium. When alignment occurs, the last 
term in (^E^ disappears, restoring the standard action for the fluid universally coupled to 
gravity. In other words, the violation of equivalence principle will be screened inside a dense 
medium, realizing an analog of the chameleon mechanism [IQ]. This picture is elaborated 
quantitatively in what follows. 

The action ( 128|) must be supplemented by the non-relativistic limit of the Einstein-aether 
action For simplicity, in the rest of this section we will restrict to the case C2 = C3 = 
C4 = and describe the results to leading order in ci <^ 1. This limit shows already the key 
features of the generic case, to be considered in Sec. 14.41 In this approximation, the action 
(|2]) at post-Newtonian order reads 

S^ = ^ d^x [40A7/' - 2^A^ + ciu'Au'] . (29) 

Additionally, we will assume Ci < Y . This is the most interesting regime from the phe- 
nomenological perspective since we do not expect to obtain bounds on the LV in DM matter 
sector that would be much stronger than the bounds on Ci. 
The equations of motion following from fl28|) and fl29|) arj^ 

(1 -Y)V' + di(f) + Y{dtu' + VWju' - V^diU^ + u^diU^) = , (30a) 

^l' = ^ , (30b) 

2M^A(p = p , (30c) 

M^ciAu' = Yp {u' - V') . (30d) 

This system is closed by adding the continuity equation for DM, 

dtp + diipV') = . (30e) 



^The case of negative Y can also be interesting from the phenomenological viewpoint if the growth of the 
instabihty induced by the negative potential is smaller or comparable to the Hubble rate. To satisfy this 
requirement the absolute value of Y must be of order Cq, cf. Eq. (|3T|) . which already sets a strong restriction 
on LV in DM. 

^°In this and only in this section dot denotes the total (material) derivative with respect to t. 
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Consider a spherical DM halo of size Rh, constant density ph and moving as a whole 
with velocity Yh with respect to the preferred frame. According to Eq. ( 130dp the aether 
perturbations inside the halo acquire the effective mass 



2 ^QIN 



Let us first assume that the halo is small, 

Rh < m;^ , (32) 

so that the range of the aether interactions exceeds the halo size. Note that this condition 
is equivalent to 

cph,s < Ci/F , (33) 

where (j)h,s ~ PhRh/^o is the gravitational potential at the surface of the halo. Then in 
f l30d|) one can neglect the M*-term on the r.h.s. and obtain for the aether field produced by 



the halo, 

u' = ^0h(x,t), (34) 

Cl 

where we have expressed the result in terms of the halo's gravitational potential 

0,(x,t) = - . (35) 

We now study the motion of a test DM particle in the field of the halo. Substituting 
into f l30a|) we obtain for the particle's acceleration 

V!p = -Y=Y+ (T^m ^ " ^^^^^'^''^^ + yiviA^H - vivid,^,] . (36) 

One observes two modifications compared to the standard case. First, the gravitational 
acceleration towards the halo is enhanced by the factor 1/(1 — F) due to the reduction of 
the inertial mass of the test particle. This effect will play a key role in what follows. Second, 
the combination in the square brackets gives rise to the velocity-dependent interaction^^! 
discussed in [13]. The latter interaction is important as compared to the change in the 
gravitational acceleration only if the halo and the particle move fast enough, 

\^H,tr>\>M. (37) 



^^Note that this interaction violates the (naive) third Newton's law: e.g. for a two-body system TO(i) V^\j + 

"^(2)^(2) 7^ ^- This is consistent with the observation [TJ that in the presence of velocity dependent forces the 
conserved momentum contains additional contributions, depending both on the velocities and the distance 
between the bodies (see also [55]). 
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We are going to see (cf. (!52|) ) that the bulk velocities of the linearized cosmological pertur- 
bations never satisfy the condition ([ST]), so that the velocity-dependent forces do not play a 
significant role in the evolution of the Large Scale Structure. However, they may influence 
the dynamics of non-linear structures such as galaxies and galaxy clusters. We leave the 
study of this issue for the future. 

Now consider the case when the size of the halo is larger than the inverse of the aether 
mass in its interior, 

Rh > . (38) 

Equivalently, 

(t)h,s > ci/Y . (39) 

Due to the Yukawa screening, the aether field is frozen at the minimum of its potential, 
= Vfl over most of the halo volume, deviating from this value only in a surface layer of 
width m~^. Neglecting this deviation, to leading order the aether field outside the halo is 

= ■ (40) 

Clearly, this is smaller than (IMj) provided ( 139|) holds. Thus the velocity- dependent interac- 
tion of a test particle with the halo is suppressed in this case. One notices the similarity with 
the chameleon mechanism ^D] that screens scalar interaction in certain scalar-tensor gravity 
models. However, it should be stressed that the motion of a test particle is still modified 
compared to GR due to the renormalization of the gravitational acceleration (cf. fl36|l ). 

Finally, instead of the test particle one can consider another large halo, such that its own 
gravitational potential satisfies f l39|l . Let us denote its velocity by Yhi- Irrespectively of the 
aether configuration outside, in the interior of this halo the aether will coincide with V^j^. 
Substituting this into (l3Uap we obtain the standard equation, 

= -^^<iyh ■ (41) 

One concludes that for large halos the deviations from GR are completely screened^. 



3.2 Jeans instability 

In this subsection we analyze how the above peculiarities of DM dynamics in the presence 
of LV affect the development of the Jeans instability. For simplicity, we will consider the 
universe filled exclusively with DM; the more realistic case including the cosmological con- 
stant, baryons and radiation will be studied in subsequent sections. Following the standard 
treatment of the Jeans instability in the Newtonian limit, we take the spherically symmetric 

""^^One may wonder if this screening mechanism can work for ordinary matter and help to relax some 
constraints on LV in the Standard Model. Unfortunately this is not the case as in the Standard Model LI is 
tested for individual elementary particles. 
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Ansatz for the background values of the fields which we denote with an over-bar. In this 
case the solution to Eqs. (15(11) reads 



H{t) x' 



Pit) 
12M2 



where the Hubble parameter H and the density p are related by 

r2_ P 



6M2 ' 



p + 3Hp = . 

These are the standard equations of GR. In particular, introducing the scale factor 



a{t) = exp 



H{t')dt' 



and integrating Eqs. 



once, we obtain the Friedmann equation 

-2 _ P 



3M2 



(42) 

(43a) 
(43b) 

(44) 
(45) 



with p oc a~^. Thus at the level of the homogeneous cosmology the model reproduceqlj GR. 
In what follows we will set the spatial curvature to zero, k = 0. 

The next step is to write down the equations for the linear perturbations. We obtain 



(1 - Y){dtV' + Hx^djV' + HV') + + Y{dtu' + Hx^dju' + Hu' 

2M^A(f) = p6 , 

M^ciAu' = Yp{-V' + u') , 

dt6 + Hx'di6 + diV' = , 







(46a) 
(46b) 
(46c) 
(46d) 



where for simplicity we denote the perturbations by the same letters as the total quantities. 
We have introduced the density contrast 6 = 6p/p. Performing the Fourier decomposition 

^ ^gikx/aw^ ^ etc. , (47) 

(48) 



restricting to the longitudinal part 



and solving fl46bl) . fl46dp for (p and v we are left with two equations 



3/72 U2 

, nzJ^,^ s + Y—{dtU + Hu) 

2 d 



[l-Y){d^6 + 2Hdt6) 



ciM^— + Yp 



a 



u 



Y pa 



dt5 



(49a) 
(49b) 



^•^As will be discussed in Sec. |4j there is actually a mismatch between the gravitational constant appearing 
in the Friedmann equation and the locally measured value (cf. ([57)) ). This mismatch is of order Ca and is 
not captured by the approximation used in this section. 
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Here we have used the Friedmann equation (H5l) to express the factor in the second term of 
fl49ap through the Hubble parameter. 

The form of Eq. ( ]49bp suggests to distinguish two regimes. For short waves, 



^ _Yp_ 

the aether perturbation reads 



2 » T75- . (50) 



Note that the inequahty f l50|) is nothing but the condition for the unscreened regime, Eq. fl32|) . 
with the size of the halo replaced by the physical wavelength a/k oi the perturbation. Upon 
substitution into fl49ap the aether gives only a small correction that can be neglected. This 
can be understood as neglecting the velocity-dependent interaction. Indeed, one estimates 
the velocity perturbation as (see (]46d|) ). 



|V| ^ « H,^^ . « ^, ,52) 

where in the third relation we used the background Friedmann equation; this is below the 
limit fl37|) . Thus for the density contrast we obtain the equation, 

d^S + 2Hd,5 - 6 = 0. (53) 

Its growing solution has the form 

'25 -Y 



6 (xf^ , 7 = - 
o 



- 1 + 



1-Y 



(54) 



Clearly, for Y > the mode grows faster than in GR where 7 = 2/3. 
In the opposite regime of long modes, 

the aether takes the form 

u = — dtS . (56) 

The corresponding contribution in fl49ap combines with the first term and restores the stan- 
dard equation for the density perturbations in GR as expected in the screened regime, 

d^6 + 2Hdt5 -^5 = 0, (57) 

whose solution is 5 oc t^/^. 
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The following picture emerges from the above results. For a mode with given k the effects 
of LV are screened as long as the inequality ( 1331) is satisfied and the mode behaves exactly as 
in GR. As the universe expands, this inequality breaks dowE0 and the mode enters into the 
unscreened regime, where its growth is accelerated up to fl54l) . Physically, this enhancement 
is due to the increase in the gravitational acceleration of the DM particles in this regime, 
see Eq. fl36|) . which in turn is due to the change in their inertial mass. We will see in Sec. [5] 
that this picture coincides with the results of the full relativistic analysis in the appropriate 
limit (subhorizon modes, matter dominated epoch). 

Before closing this section, we mention that a straightforward analysis of the transverse 
(vector) modes of the perturbations shows that these modes decay in the same way as in GR 
both in the long- and short-wave regimes. We will not consider the vector modes in what 
follows. 



4 Lorentz violating dark matter: relativistic fluid 

To systematically develop the consequences of the aether - DM interaction beyond the 
Newtonian limit it is convenient to use a relativistic fluid description of DM. The most 
suitable for our purposes is the pull-back formalism, which enables to describe perfect fluids 
with arbitrary equation of state using an effective action, see [H] for a review and early 
references. Since the reader may not be familiar with this formalism, we give a succinct 
review of it in the next subsection, closely following the presentation of [29j. We will specialize 
in the case of cold DM with vanishing pressure and present an effective action describing its 
interaction with the aether. In Sec. 14.21 we use this action to derive the equations of motion. 
Then in Sees. 14.31 14.41 we analyze the consequences of the resulting model for the background 
cosmology and derive the equations for linear perturbations around Friedmann-Lemaitre- 
Robertson-Walker (FLRW) solutions. 



4.1 Effective action for Lorentz violating fluids 

In the pull-back formalism, the fluid elements are labeled by three scalar fields (p^{x), I = 
1, 2, 3. An essential property of a perfect fluid is the invariance of its dynamics when its 
elements are moved around without changing their volume. This means that the description 
in terms of the fields ^p\x) must be invariant under volume-preserving reparameterizations 

^ ^\^') , det 1^ = 1 . (58) 

The scalar object with the minimal number of derivatives invariant under (1351) that can be 
constructed from ip^ and the metric Qf^u is 

B = -det B^'^ , (59) 
^^Recall that p decays as a^^ so that the r.h.s. of ([55)) decreases faster than the l.h.s. 
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where 

B'' = g^'dyd.if' . (60) 
The Lagrangian for the perfect fluid is thus an arbitrary function of this quantity, 



S fluid = - I d^xy^fiB). (61) 
To make contact with the standard quantities of hydrodynamics let us consider the vecto'^ 



d.^'d^^'dp^^'eijK . (62) 



This vector has unit norm, is future directed and has the property that the fluid labels (p^ 
remain constant along its integral curves 

yt^Q^^i = J = 1,2,3. (63) 

Thus it is identified with the fluid velocity. It follows that independently of the equations of 
motion 

V^{^v^) = ^, (64) 



which implies that v B should be interpreted as the conserved number density n of the fluid 
elements. Other useful relations are 



BijdpLp^dyip^ = Qf^u - Vf,Vu , dyip^dpLp-^ d^Lp^eijK = -\/-gB v^t^upa , (65) 

where Bjj is the inverse matrix of B^"^ . The energy-momentum tensor (EMT) following 
from 0611) can be cast into the standard hydrodynamic form 

Tpu = (p + p)vpV^ - pg^y , 

where the density and pressure are given by 

p = /, p = 2f'B~f. (66) 

It is conserved due to the equations of motion for the fields ip^ , 

V^(f(5)55,jVV) = 0. (67) 

For cold DM p = implying that / takes the form, 

fdm = my/B = mn , (68) 

where m is the mass of the fluid elements. 



^'^The antisymmetric pseudotensors are defined as e^i^a _ —60123 = — 1, £123 — 1- 
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It is clear how to generalize the action ( 16T|) to include the effects of LV. In the presence 



of the aether the action can also depend on the invariant scalar product u^v^, 



S 



fluid 



(69) 



Once the microscopic theory of DM is specified the function / can in principle be computed. 
In this work we adopt a phenomenological approach and consider a generic function satisfying 
reasonable conditions. A natural assumption for DM is that the Lagrangian is proportional 
to the particle number, 

fdmiB, = , (70) 

where the function F is subject to the normalization condition -F(l) = 1. It is easy to 
convince oneself that most cold DM scenarios are described by functions of the form (I7CT]) . 
In particular, this assumption is valid for DM composed of well-separated particles that 
weakly interact with each other. The function F then coincides with that appearing in the 
one-particle action considered in Sec. |3l 



4.2 Equations of motion 

The equations of motion for the DM fluid are obtained by varying fl69|) with respect to ip^ 
and read 



BBjjV^<f'{F - {uxv^)F') + F'-^u,dx^-'dpip''eijK 



. 



(71) 



A more conventional form involving only hydrodynamic variables emerges after contracting 
with da(p^ and using Eqs. f l65|) . 



-V,, 



P[dm]{{F - {uy)F')v^v„ + F'v^u^y^ + p^d^F'v^^V.u^ = 



where we introduced 



Pldn 



ni\fB . 



(72) 



(73) 



The equations of motion for the aether are slightly different depending on whether one 
considers a generic aether or a khronon for the LV sector. In the former case, from the 
variation of the combination of (12]) and (ESI) one finds 



Vf,K^^ - c^aPVyUp - lu 
where we made use of the following notations. 



_ Pldm] , _ 



(74) 



aa = U VxUn ■ 



(75) 
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Contracting (17^ with and using the constraint ([T]) we obtain the expression for the 
Lagrange multipher, 

^-^{u^v^F'. (76) 



When substituted back into (1741) . this yields, 



(77) 



where 



(78) 

is the projector on the directions orthogonal to u^. In the khronon case the constraint is 
identically satisfied and thus the Lagrange multiplier is absent. Varying with respect to cr, 
we find the equation of motion which is essentially the divergence of fl77p . 



. 



(79) 



The EMT of DM is obtained by varying fIBI?]) with respect to the metric. In the case of 
the Einstein-aether we obtain 



T[dm]fiu = P[dm] {F - {uxV^)F')v^Vy . 



(80) 



In the khronometric model there is an extra contribution coming from the variation of 
when the metric is varied, 



Tr 



[dm] iiu P[dm] 



{F - {uxv^)F')v^v, + {uxv^)F'u^u, 



The difference between the previous expressions disappears once we consider the total EMT 
of the aether / khronon — DM system. This reads in both caseff], 



2V,K\^u,) - Vx{K^,u)y') - V.{K\^u,)) + V.(ir(/n.)) 



+ P[dm] ({F- {uxv^)F')v^v, + {uxv^)F'u^u, 



12) 



where in the aether case we have substituted the Lagrange multiplier from (I76|) . The first 
three lines in (18 2 p correspond to the pure aether EMT, while the last line represents the 
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The indices in round brackets are symmetrized as K(^^^) = ^{Kf^i, + K^f,) , etc. 
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contribution of DM and DM — aether interaction. Note that LV leads to the departure of 
the latter contribution from the perfect fluid form. 

For the rest of the Universe we will consider the standard scenario with baryons, radiation 
and a cosmological constant^. Those possess EMT's, 

T[b]fiu = P[b]V[b]^iV[b]u , Tl^]^lu = P[7] ^g^HM^W^^ ~ ' ^[aIm^^ = P[A]9^lu , (83) 

where p[s], V[s]^ is the density and 4- velocity of the corresponding component. The equations 
of motion for these components are given by the covariant conservation of the EMT's fl83|) . 
We neglect the energy exchange between baryons and radiation, as well as the effects due to 
neutrinos. 



4.3 Background cosmology 

The spatially homogeneous and isotropic Ansatz reads 

ds^ = a^(r)(dr^ - dx^) , uq = V[s]o = a{T) , m = vis]i = , p[s] = P[s]ir) , (84) 

where we have introduced the conformal time r. It is straightforward to check that with 
this Ansatz Eq. ( I77p for the aether, or alternatively Eq. (17^ for the khronoij^. is identically 
satisfied. Substituting f lH^ into fl72|) we find 

Pldm] oc , (85) 

so the DM density behaves in the same way as in the standard FLRW universe. Using (152]) 
we obtain the Friedmann equation, 

-I = o [P[dm] + P[b] + P[7] + P[A]) , (86) 

Qi o 

where dot denotes henceforth differentiation with respect to the conformal time. This has 
the same form as in GR, with the gravitational constant renormalized due to the contribution 
of the aether (khronon) EMT in the Einstein's equations (cf. [32]). 



(87) 



Note that Gcosm differs from the Newton constant ([9]) measured in local experiments. This 
difference is constrained by the BBN considerations [32] 

\Gcosm/GN — 1| < 0.1 , 

which places a rather mild bound on the aether parameters. 



we want to isolate the effect of LV in the DM sector, we will not consider other possibilities for dark 

energy. 

^^The khronon field itself can be taken as an arbitrary monotonic function of time, e.g. a = t. 
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4.4 Linear perturbations around FLRW 

We now turn to the analysis of the hnearized cosmological perturbations. We focus on scalar 
perturbation^^ and choose to work in the conformal Newtonian gauge, 

ds^ = a(r)2 [(1 + 20)dr2 - 6ij{l - 2^)dx'dx^] . (89) 

Since the scalar sectors of the Einstein-aether and khronometric theories are equivalent at 
the level of perturbations, it is enough to study the khronometric case with the parameters 
dH]). For the perturbations of the aether vector, the velocities of various components and 
their densities we write 



Ui 



a diX , V[s] i = a diV[s] , uq = V[s] o = a (1 + 



P[s] = Pis] + Sp[s] 



(90) 



where [s] = {[dm], [b], [7]} and the overbar denotes the background values. Substituting 
these expansions into flS^ and (1551) we obtain the set of linearized Einstein's equation^. 



2Aip - m{2 + aB)i) - aA0 + aAx + aH(l - i3)Ax 



^(5P[dm] + 5p[b] + 5p[^] + 2(j)[p^d'm] + P[b] + P[7] + P[A]) = , (91a) 



'0 

.2 



2 + aB Ml{2 + aB) 



P[dm]{Yx + (1 - Y)V[dm]) + P[b]V[b] + gPbFH 



+ 7/(0 + 2^) + (27^ + 7/2)0 + 

0-^-/3(x + 2-Hx) = O, 
where we defined 



2A(0 - ^) + aBAjx + 2Hx) a^^p^^j 

3(2 + aB) 3(2 + aB)M^ 



B = 



(3 + 3X 



(3 + X 



a 



a 



a 



Y = F'(1] 



(91b) 



(91c) 
(91d) 



(92) 



^^As we mentioned in Sec. 13.21 the vector perturbations decay as usual inside the horizon, so we do not 
expect any significant effect due to thcni. For the analysis of effects at cosmological distances in the Einstein- 
aether model with Lorentz invariant DM see [32]. Tensor modes behave as in GR with a modified speed of 
propagation [T^ [T5] . 

^"In these and all subsequent equations dot denotes partial derivative with respect to r. 
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We will assume in what follows that the parameters B and C are both of order one. The 
linearized khronon and DM equations are 



'0 



X 



aM^ ^M- 0-^(1 + i3)0-i37/^ = O, (93a) 

Sp[dm] + 3'H5p[dm] - P[dm]{^V[dm] + 3^) = , (93b) 
V[dm] + 'HV[dm] + ^ (X + '^X) - " ^ • 

From f l93ap we see that the parameter C has the physical meaning of the velocity of the 
khronon waves. In what follows we will often use the integrated form of the last equation, 

''^'^'-^ = -^^"^ a{r){l-Y) j dr'a(r')0(r',x). (94) 

If the gravitational potentials are known, this can be used to eliminate v^dm] from fl93al) and 
obtain an equation for x- We observe that, as in the Newtonian limit, LV in DM is governed 
at the linear level by the single parameter Y . Following the discussion in Sec. [3l we will 
assume it to lie in the range < F < 1. 

For baryons and radiation we have the standard hydrodynamic equations, 

+ 3'H5p[;,] - p[b] (At;[5] + 3^/') = , + Tivy^ -0 = 0, (95a) 

5p[^] + m5p^^^ - ^p[^] [Avy^^ + 3^^) = , -^-0 = 0. (95b) 

The equations presented above completely determine the evolution of the linearized cosmo- 
logical perturbations. Their exact solution can be obtained only numerically. In the following 
section we find the approximate analytic form of the solution in various dynamical regimes. 
This will enable us to determine the qualitative influence of LV on the power spectra. The 
numerical analysis is postponed till Sec. Ei 



5 Cosmological perturbations: qualitative analysis 

Consider the mode with a given conformal wavenumber k. As usual, its evolution will be 
different depending on whether k is smaller or larger than the expansion rate 1-L. Coupling 
between DM and khronon gives rise to a second dynamical scale crucial for the mode evolu- 
tion. Substituting Eq. ( 194|) into Eq. fl93ap . we notice that the dynamics of the mode depends 
on the ratio between k and 

2 \ 1/2 
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This new scale is directly related to the density of DM and determines the critical wavenum- 
ber below which the effects of LV are screened (see discussion in Sec. [SD- Factoring out the 
explicit time dependence we can write, 



where Hq, Qdm are the present-day Hubble constant and DM density fraction; the current 
value of the scale factor has been normalized to one, a(ro) = 1. Note that under the 
assumption a '^Y the present screening scale ky^ is parametrically higher than the Hubble 
rate Hq. The hierarchy between ky and "H persists during the whole matter dominated epoch 
and most of the radiation domination, see Fig. [TJ Thus, the time evolution of a mode with 



Figure 1: Time dependence of the dynamical scales determining the evolution of the cos- 
mological perturbations. Tgg is the time of radiation-matter equality and tq — the present 
time. See the main text for other notations. 

given k can be divided into the following regimes: 

(a) super horizon modes, k kH, radiation dominated universe, 
(a') superhorizon modes, k kH, matter dominated universe, 

(b) subhorizon modes with % <k < ky, radiation dominated universe, 
(h') subhorizon modes with < k < ky, matter dominated universe, 

(c) subhorizon modes with < ky < k, radiation dominated universe, 
(c') subhorizon modes with H < ky < k, matter dominated universe. 




(97) 




T, conformal time 
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To find the solution of the linear equations we use the following strategy. We first consider 
the formal limit a, /3, A — )• suggested by the constraints on the khronon (aether) parameters 
listed in Sec. [H at the same time we keep the parameter Y fixed. From fl93ap one sees that 
in this limit the khronon is tightly coupled to the DM velocity, x = V[dm]- Substituting 
this into the rest of equations we find that they take the standard GR form. Thus, to this 
approximation all perturbations, in particular the gravitational potentials and ip, behave 
as in GR. At the next step we allow for finite a, P, A (assuming all of them to be of the same 
order) and substitute the GR values for 0, ip into Eqs. (jHSD to find the khronon and DM 
fluctuations. Finally, we find the corrections to 0, ip by inserting x, V[dm], Sp[dm] into the 
Einstein's equations (19T|) . 

This logic has the caveat that in the limit a, /3, A — )■ the screening scale ky diverges. 
Thus, the above perturbative scheme is not applicable in cases (c) and (c') where the mo- 
mentum k is parametrically large. Those cases will require a different approach. 



5.1 Regime (a): k <7i during radiation domination 

For the superhorizon modes, we can neglect all terms with spatial Laplacians and consider 
to leading order = 0^ = const. At the radiation domination stage the scale factor depends 
linearly on the conformal time. 



a r 



(98) 



where Vt^ is the present fraction of the radiation density. Substituting this into Eqs. fl93ap . 

and setting the constant of integration in (IMI) to zero, which amounts to the absence 
of initial velocitied^, we obtain, 



X + — + 
r 



\ ii^ 



X 



'l + B) k^, 
T 2A. 



7 J 



^7 • 



This equation has a solution 

From 0941) we find that it corresponds to 



(99) 



(100) 



V[dn 



X 



-T . 



(101) 



This is the standard expression for the perturbation of DM velocity in the adiabatic mode 
and we see that the khronon evolves coherently with the DM fluid. Note that this solu- 
tion is completely insensitive to the parameter Y setting the strength of LV in DM. The 
corresponding DM density contrast in the adiabatic mode is [13] (see also Appendix Rl) . 



Si 



[dm] 



Sp[dm] 
P[dm] 



(102) 



'^Even if the initial velocities are present, their effect decays with time. 
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Inserting f llOip into the Einstein's equations (pTj) one can find the corrections to the gravi- 
tational potentials induced by the khronon. Clearly, these corrections are suppressed by the 
parameters a, (3, A and we are not going to analyze them in detail. Let us just note that 
Eq. fl91d|l implies the appearance of a small anisotropic stress, 

i^.f. ,103) 

To determine whether the adiabatic mode fllOll) is an attractor, we search for the solution 
of the homogeneous equation, i.e. Eq. (199]) with vanishing r.h.s. This can be found in terms 
of Bessel functions. 



where Xi^l arbitrary coefficients. Recall that we restrict to < F < 1, in which case 

kYfi is real and thus xioL oscillates and decays at late times as r"'^/'*. In this case the 
adiabatic solution f llOOp is an attractor. For negative Y, ky^ is imaginary and the solution 
f ll04p contains an exponentially growing part signaling the instability of the system. We 
have already encountered this instability when we studied the Newtonian limit in Sec. |3l 
For the rest of the paper, we will only consider positive values of Y . 

At early times, r <^ A^/kyQ-, the first term in fll04p evolves as r"^ with 



q = . (105) 

For i3 < — 1 it grows faster than the adiabatic mode fllOOp and, depending on the initial 
conditions, may dominate over it. This corresponds to an intrinsic instability of the khronon 
(not to confuse with the instability discussed in the previous paragraph) [121 [22] • It is 
curious that the coupling to DM cuts off this instability at r ~ A^/kyo, and makes the 
apparently growing mode decay. In this paper we will assume that the initial conditions for 
the cosmological perturbations are sufficiently close to adiabatic, so that this mode never 
becomes dominant. 

5.2 Regime (a'): k <% during matter domination 

During matter domination the scale factor grows as 

a(r) = A^T^ , = (fi^^ + n,)Hl/A , (106) 

where VL^ is the baryon density fraction today. At zeroth order, the superhorizon amplitude 
of the gravitational potential is now 0^ = ^07- After substituting this into Eqs. f l93ap . flM|) 
and neglecting the terms with spatial Laplacians one obtains the following equation: 



X + — + 

T 



k"^ H 

2(l + 3i3) + ^ 



r2 



2{l + B) + 



k 



2 1 

Yfl 



3Ar. 



— . (107) 

T 
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Here we have neglected the integration constant in ( IMj) which corresponds to a decaying 
contribution. The adiabatic mode 



V[dm] = X = (108) 
is a solution. To check its stability consider the solution of the homogeneous equation, 



{a') I N (a') r+ , (a!) r_ 



^ 2 



3 ± W 1 - 24i3 ^'^ 



An 



(109) 



a') 



For all parameters of interest the expression under the square root is negative and 
describes decaying oscillations. The adiabatic mode f llOSp gives rise to the anisotropic stress 

^ = M . (110) 



Comparing to (11031) and taking into account the ratio between 0^ and 0^ we see that the 
absolute difference (0 — ■?/') actually stays constant at the transition from radiation to matter 
domination. 



5.3 Regimes (h)^ (c): k > Ti during radiation domination 

The gravitational potentials decay rapidly inside the horizon during the radiation domina- 
tion, so we can neglect them in Eqs. fl93al) . flM|) . The latter then takes the form. 



Y 



\drn\ 



x + 



1-Y^ r{l-Y)' 



:nv. 



where v^^^ is a constant. By matching fillip approximately to fllOip at horizon crossing 
fcr ~ 1 we find ^ ^^/2k'^. Substitution into flMaj) gives 



X + — + 

T 



1.2 1 



X 



;il2) 



where we have neglected the terms related to the Hubble parameter T-L in the square bracket. 
For k ^ kyfi/ \fAlfT (regime (hj) the solution reads. 



yi'^ 



2k 



Y,0 



(b) 
X2 



2k 



Y,0 



(113) 



This is a sum of a monotonically decreasing mode and damped oscillations. The same 
structure is inherited by V[dm] through Eq. fillip and eventually, due to Eq. f]93bp . leads to 
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the standard logarithmic growth of the DM density contrast (the LV effects are screened), 
accompanied by damped oscillations, 

3(f) (j) 

5[dm] ~ — ^ — ^ log {kr) + (damped oscillations) . (114) 



For k ^ ^y,o/ \/ ^■y'T (regime (c)) the formal perturbative expansion in a,/3,A breaks 
down, as the corrections proportional to ak"^ in the Einstein's equations (19T|) may become 
large in this regime (recall that kY,o is inversely proportional to y/a). Nevertheless, we can 
rely on the fact that during the radiation domination the DM perturbations are subdominant. 
It is reasonable to assume that this is also true for the khronon fluctuations which are coupled 
to DM. Then the dynamics are still dominated by the sound waves in the hot plasma that 
wash out the gravitational potentials. Therefore, Eq. f illip for the velocity of DM and 
Eq. flll2[) for the khronon still hold. Moreover, in the latter equation one can neglect the 
second term in the square brackets and the r.h.s. and obtain. 



(c) sm [LkT) (c) cos [LkT) 
x(^) = Xi + X2 • (115) 



sin {Ckr) (c) cos {Ckr) 

T T 

Using f lllll) . (]93bp we obtain for the density contrast. 



= const - ^(Y^ + ^^^Si(CA;r) + ^M!ci(Cfcr) (116) 

where Si(x) and Ci(x) are the integral sine and cosine, 

sin^ , cos^ — 1 



r sin f r'-' 
Si(x) = / di^, Ci(x) = 7 + ln(x) + / 
Jo s Jo 



Notice that the coefficient in front of the logarithmic term has been renormalized by the LV 
effects. The last two terms again describe damped oscillations. 

The oscillations in the DM density and velocity field found above could potentially provide 
an interesting signature of the model. Unfortunately, the numerical solution in Sec. E] shows 
that they are rather weak and can leave a noticeable imprint only at very short wavelengths 
that are presently in the non-linear regime. It is unclear whether such effects can be observed. 



5.4 Regimes (b'), (c'): k > Ti during matter domination 

We finally consider the evolution of subhorizon modes at the matter dominated epoch. If 
one neglects the effects of LV, the Newtonian potential stays constant, (j) = 4'm,ki where the 
subscript k indicates the dependence on the wavenumber. Neglecting in flMjl the integration 
constant that gives a rapidly decaying contribution and substituting the resulting expression 
for the velocity into fl93ap we obtain 



X+ — + 
r 



X 



2{l + B) + 



3 Am 



T 



(117) 
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For k <^ kyfl I V Am T (regime {h')) the solution is 



X 



6/lm(l + i3) + A;: 



2 



2 



'rin,lc ' I Al ' 



+ , p± = --±W---±- . (118 




The definitions of kyfl, Am in (197|) . (11061) imply that fcyoMm ~ is parametrically large. 
Thus the last two terms in the above solution describe quickly decaying oscillations and one 
is left with the attractor behavior (cf. (llOSp ). 



V[dm] ~ X ~ ■ (119) 

The velocity perturbation has exactly the same form as in GR, so we again conclude that 
in this regime the LV effects are essentially screened. The corresponding density contrast is 
found from fl93bl) and exhibits the standard quadratic growth in time, 



S[dm] ~ -^-^<Pm,k . (120) 

It is instructive to go one step further and consider the corrections to the gravitational 
potentials induced by the khronoiJ^. From Eqs. (I91cp . f lQldp we find 



^//i,\\/i.\2\ ; J /i ^/^ I ^ ro I \\n \2 



i + ^W + x){kTy , ^ = ^^ji-^ + -(i3 + x){kTy] . (121) 



These expressions are valid up to linear order in A. We observe that the corrections, 
though small, grow with time. Eventually, when ky red-shifts down to k, they become of 
order Y: this follows from the estimate (fcyr)^ ~ Y/a. The suppression of the corrections by 
the small parameters A disappears, which suggests the break down of our perturbative 
scheme by the end of the regime (b'). Note however, that the anisotropic stress remains 
constant and given by Eq. flllOp during this stage. 

For k ^ ky^Q/^/A^T (regime (c')) one can show that the corrections to the gravitational 
potentials are not proportional to the expansion parameters a, (3, A, so the perturbative 
calculation fails. Thus we have to solve the coupled Einstein-khronon-DM equations self- 
consistently. One notices that fl91cp . flQldl) . fl93al) . fl93cl) form a closed system of equations 
for 0, ijj, X, V[dm]- The educated guess of a power-law behavior yields the growing mode, 

^ = = ^m,k , (122a) 



^^There are also corrections to S\^cim]i 'V[dm] the same order which we do not consider here. 
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where 

5 /25 QY 



and we have neglected the corrections of order a, P,X. The mode normahzation (f)m,k can be 
determined by matching fll22p to the solution before crossing the screening horizon k = ky- 
Substituting fll22cp into (]93b|) we find the behavior of the DM density contrast, 



(- k (f)m,k /'10/l^ 

= -(i-r)(x + 3)(x + 2) " • 

One observes that 6idm] grows faster than in GR, implying the accelerated growth of structure. 
The power-law ( 11241) coincides with the results in Sec. 13.21 (see Eq. f lMl) ) if we neglect the 
baryon contribution, i.e. set = 0, and take into account the relation between the conformal 
and physical time, oc t. 

It is interesting to study the perturbations of the baryonic component. From fl95al) we 
find, 

''^'= p;;r = "(-+3)(x+2)" • ^'''^ 

Thus the baryon density contrast grows at the same rate as that of DM, but with an overall 
amplitude suppressed by a y-dependent factor. 



S[dm] 



1-Y. (126) 



This behavior has a simple physical explanation: baryons do not feel the enhanced grav- 
itational acceleration experienced by DM. Therefore it takes them longer to react to the 
inhomogeneities of the gravitational potential and their density perturbations stay behind 
those of DM. 

The anisotropic stress is found from fl91dp and reads. 



"■^ P-^J^- (127) 



' C^A^kH"^ 

Notice that the second factor on the r.h.s. is always smaller than 1 in the considered regime, 
so that the anisotropic stress is still suppressed by /3. Moreover, it decays with time as r~^. 

Finally, it is straightforward to check that the solution f ll22p . f ll24p . fll25p satisfies the 
remaining Einstein's equations fl91ap . (191bp . 



5.5 Qualitative analysis of the power spectra 

The previous analysis allows us to sketch the qualitative picture of the effects due to LV in 
DM in the observed power spectra of DM and baryon density perturbations. From Fig. [1] we 
see that the modes with different k go through a different sequence of the dynamical regimes 
described before. 
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log k 



Figure 2: Schematic representation of the matter power spectrum in the LV dark matter 
model (upper curve) compared to the power spectrum in ACDM (lower curve). The same 
qualitative behavior is common to the dark matter and baryons. The quantities ky^, ky^eq, 
X are defined in fp7|) . fll28p . (11231) . The figure corresponds to the case ky^Q < l/Teq- 



• The modes with k < ky^ always stay in the regime where the LV effects are essentially 
screened. In this range the power spectrum has the same dependence on k as in ACDM. 

• The modes with ky^Q < k < ky^^g, where 

kY,eq = ^. ^ = ky^o\ 7^ , (128) 

is the screening scale at the radiation-matter equality, enter into the regime of un- 
screened LV at the matter dominated stage. Let us denote the moment of the "screen- 
ing horizon" crossing by Ty{k), 



Ty{k) = 




Between ry and tq the density contrasts exhibit the anomalous growth fll24p . (I125p . 
Thus they are enhanced by a factor (ro/ry)^ leading to the increase of the power 
spectra by a factor oc fc^^ in this range of momenta. 

• Finally, the modes with k > ky^eq enter in the unscreened regime already during radia- 
tion domination. Modulo the small damped oscillations mentioned in Sec. 15. 3[ the net 
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effect of LV on tliese modes is tlie overall enhancement by the factor (To/rgg)^. This 
factor does not depend on the mode momentum, so the slope of the spectrum remains 
as in ACDM. 

This picture is summarized in Fig. |2j For the sake of the schematic representation we have 
assumed a flat spectrum for the initial perturbations, neglected the logarithmic growth of the 
perturbations during the radiation domination together with the effects of the cosmological 
constant. Note that the change in the slope of the power spectrum depends only on the 
parameter Y describing LV in DM, while the range of scales where this change occurs is 
determined both by Y and the khronon parameters a, /3, A. 



^ 1-Y 
CO 

to 



log k 



Figure 3: Ratio between the density contrasts of baryons and dark matter (qualitative plot). 

Figure [2] corresponds to the situation when the ratio Y/ a is moderately large, so that 
kvfi is smaller than l/rgg, the Hubble rate at the radiation-matter equality that sets the 
position of the maximum of the power spectrum in ACDM. In this case the change in the 
slope of the power spectrum leads also to the shift of its maximum. Indeed, in the log-log 
scale the standard ACDM power spectrum in the vicinity of the maximum can be written as 

logP^^^^^ = Ai-A2(logA;re,f , 

where 2 are some constants. The additional enhancement due to LV adds a linear contri- 
bution on top of it, 

logP = logP^^^^^ + 2xlogfc . 
Differentiating this expression with respect to k we obtain the new position of the maximum. 



log kmax = log 1/ Teq + >i j A2 

We conclude that the shift is linearly proportional to k. 



(129) 
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If kyfi > l/Teq, which happens for an extreme hierarchy between Y and a, the position 
of the maximum of the power spectrum is not modified. However, the change in the slope is 
still present providing the signature of LV. 

The second qualitative effect of LV in DM is the appearance of the baryon — DM bias 
(11261) for the modes undergoing the anomalous growth, i.e. all modes with wavenumbers 
larger than kyfl- For smaller k the density contrasts of the baryons and DM are equal. 
This results in the scale dependence of the bias shown in Fig. [31 It would be interesting to 
understand how this effect could be constrained by observations. 

Finally, at non-zero value of the parameter /3 the model is also characterized by the 
anisotropic stress. Its present-day value is independent of k for k < ky^ where it is given by 
Eq. (11 101) . At larger wavenumbers it falls off as according to Eq. (I127p . 

The qualitative analysis of this section is confirmed by the numerical calculations, to 
which we proceed now. 

6 Cosmological perturbations: numerics 

The study of the effects of LV in the DM sector yielded the modified linearized equations 
for cosmological perturbations (PTi) . (jnSD, dHSD- In this section we will solve these equations 
numerically for different sets of LV parameters and confirm the modifications in the growth 
of perturbations uncovered in Sec. [51 The density fractions of the cosmological constant, 
DM, baryons and radiation are taken to be 

= 0.75, ndrr. = 0.2, fib = 0.05, n^ = 5- 10"^ 

We neglect the interaction between baryons and radiation as well as the effects of neutrinos. 
As our goal is to compare the evolution of perturbations in models of LV DM with ACDM, 
we also find the evolution of perturbations in ACDM within the same approximation^. The 
details of the numerical procedure are presented in Appendix [Al 
The first set of parameters we consider is 

a = 0.02 , /3 = 0.01 , A = 0.01 , F = 0.2 . (130) 

This choice satisfies the PPInI^. BBN and gravitational radiation constraints [TH [18]. The 
initial conditions correspond to the adiabatic mode, and for illustration purposes the initial 
value of is normalized to 1 for all modes. Figure [H shows the dependence of the Newton po- 
tential on conformal time for the choice (I130p and for ACDM for several modes with different 
values of momentum. The difference (0 — tp) is also shown. We see that perturbations of the 
gravitational potential are enhanced at late times as compared to ACDM. The enhancement 
is stronger for shorter modes that enter earlier into the regime where the LV effects are not 
screened. Note that the present-day screening scale corresponding to the choice (I130p is 

^■^The corresponding equations follow from (IMI) . (j93bl) . (|93cl) by setting Y = a = l3 = X = Q,x = V[dm]- 
^■^Recall that a = 2/3 case avoids the PPN bounds. 
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Figure 4: Dependence of the gravitational potentials on conformal time for modes with 
momenta k/Ho = 2, 4, 25, 300. The LV parameters of the model are a = 0.02, /3 = 0.01, 
A = 0.01, Y = 0.2. The standard ACDM result is shown for comparison. Dashed vertical 
lines mark the present conformal time tq = 3.5Hq^. 



kyfl = 2.77 Hq. On the other hand, the difference between the two gravitational potentials 
{(p — ip) which is initially of order 10~^ decreases once the mode enters inside the screening 
horizon. Note that the overall amplitude 10~^ for long wavelength modes agrees with the 
estimate {(p — ip) / (j) ^ (3 obtained in the previous section. 

In Fig. O we plot the total (DM plus baryons) density contrast in the present model 
divided by its value in ACDM. As expected from the analysis of Sec. [5] (or simply from the 
enhancement of the gravitational perturbations), the growth of structures increases at recent 
times in a wide range of momenta. The relative effect is stronger for shorter modes and can 
be large for our choice of parameters. 

In Sec. 15.31 we found that the subhorizon DM perturbations exhibit damped oscillations 
during the radiation-dominated epoch. To study this effect we plot in Fig. E] the early-time 
behavior of the longitudinal component of the DM velocity (see Eq. f l90|) for the definition) 
and the DM density contrast; the latter is normalized to the density contrast in ACDM. 
One observes that the oscillations are rather weak in the velocity and are almost completely 
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Figure 5: Ratio of total density contrasts between the model with LV DM and ACDM versus 
time for several values of the mode momentum k. The parameters are the same as for Fig. |H 
To is the present conformal time. 




Figure 6: Time-dependence of the longitudinal component of the DM velocity (left panel) 
and the DM density contrast (right panel) for several subhorizon modes at the radiation- 
dominated epoch. The dashed vertical line marks the time of radiation - matter equality, 
Teq ~ 0.023i/g"^. The parameters of the model are the same as for Fig. |H 



washed out in the density contrast. Besides, for the velocity they are perceptible only in 
very short modes that have by now entered into the non-linear regime. Thus it is unclear 
whether the oscillations can have an impact on any observable. 

We now consider the effects of different sets of parameters on the spectra of perturbations 
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Figure 7: Matter power spectrum for several values of the parameters listed in Table [H The 
case of ACDM is shown for comparison. 

at the present moment of time. The different choices are hsted in Table [1] together with the 
corresponding values of the screening scales fcy^o and ky^eg (see Eqs. (P7|) . f ll28p for definitions). 
All parameter choices are consistent with the gravitational tests described in Sec. HI The 
initial spectrum is taken to be flat with the same normalization in all cases. 





a 


(3 


A 


Y 


kyfl (h Mpc-^) 


kY,eg (h Mpc-1) 


a 


2 ■ 10-2 


10-2 


10-2 


0.2 


9.2 ■ 10-4 


6.5-10-2 


b 


2 • 10-^ 


10-4 


10-4 


0.2 


9.1 ■ 10-3 


0.65 


c 


2 • 10"^ 


10-4 


10-4 


0.02 


2.6 ■ 10-3 


0.18 


d 


10-^ 





10-^ 


0.2 


0.41 
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Table 1: The values of the parameters used in numerical simulations. 

The comparison between the matter power spectrum in the LV models and in ACDM 
is shown in Fig. [71 The left panel shows the cases when the present screening momentum 
kyfl is lower than k^ax — the position of the power spectrum maximum. We clearly see the 
change in the slope of the spectrum in the interval ky^Q < k < ky^eq accompanied by the shift 
of the position of the maximum. The effect is significant for values of the parameter Y as low 
as a few per cent, which suggests that these values can be tested observationally. The right 
panel shows the situation when fcy o is larger than femax; corresponding to very small values 
of the khronon parameters a, /3, A and relatively strong LV in DM, see Tabled] The position 
of the maximum does not move in this case but the change in the slope is still visible. 

Figure |8] shows the ratio between the amplitudes of perturbations in the baryonic and 
DM components. As expected from the analytic considerations of Sec. |5l this ratio drops 
from 1 at A; < ky^ to (1 — Y) at larger momenta implying a scale dependent bias between 
baryons and DM. 
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Figure 8: Ratio between the amplitudes of perturbations in the baryonic and dark matter 
components. The curves correspond to the parameters hsted in Table [TJ 



0.1 




Figure 9: Relative difference between the two scalar gravitational potentials for several 
choices of parameters listed in Table [TJ 



Finally, Fig. |9] presents the /c-dependence of the "relative anisotropic stress" — the dif- 
ference between the two gravitational potentials and ip in the conformal Newton gauge, 
normalized to 0. We observe that at small momenta it has a plateau with the magnitude 
set by the parameter (5. The plateau extends up to /c kyfii beyond which the anisotropic 
stress drops as the approximate power-law All this is in agreement with the analytic 
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estimates of Sec. [51 Note that the persistence of the anisotropic stress up to relatively large 
momentum, kyfl Hq, is a peculiar signature of the present model that contrasts with the 
more common situation where the anisotropic stress quickly decays for subhorizon modes. 
Besides, the anisotropic stress is present at early times even for adiabatic initial conditions. 
This distinguishes the present model from Lorentz-invariant models of modified gravity, such 
as f{R) IB] or DGP |15] theories, where the anisotropic stress is generated only at a recent 
epoch [in]. 

7 Summary and discussion 

In this paper we have studied the possibility to test the Lorentz invariance of dark mat- 
ter (DM) with cosmological observations. Our description is based on the Einstein-aether/ 
khronometric model, which provides an effective description of Lorentz invariance violation 
(LV) in gravitational relativistic theories. In those models, LV is encoded in a new field 
(aether) whose expectation value determines a local preferred frame. We considered DM 
as described by a pressureless fluid, and generalized its dynamics to include the LV effects. 
Those effects amount to different couplings between DM and the aether field. For cosmologi- 
cal perturbations in the linear regime all the LV effects in the DM sector can be summarized 
in a single parameter Y. This constant is to be added to the parameters of the aether sector 
(see Eqs. ([2]), ([3]) and ([8])), which are constrained by local tests of gravity. We considered 
the Newtonian limit of the model and demonstrated that LV implies modification of the 
inertial mass for small DM halos thus leading to the violation of the equivalence principle. 
For large halos this effect is screened by a variant of the chameleon mechanism. Addition- 
ally, we pointed out the presence of a velocity-dependent interaction; however, its role in the 
evolution of the Universe at large scales is negligible. 

The homogeneous expansion history of the Universe for LV DM was found to be exactly 
the same as in ACDM. However, the evolution of linear cosmological perturbations presents 
three major effects permitting us to distinguish between the two scenarios. The first effect 
is the accelerated growth of inhomogeneities for the modes affected by the violation of the 
equivalence principle. These are the modes that are short enough so that the chameleon 
mechanism does not switch on. This effect eventually leads to the increase in the slope of 
the matter power spectrum with respect to ACDM in a range of momenta. The enhancement 
depends only on the parameter Y, whereas the range of momenta is determined also by the 
aether parameters (cf. Fig. [7j). The second effect is the appearance of a new bias between 
the fluctuations of dark and ordinary baryonic matter. Importantly, the bias exhibits scale- 
dependence already at the linear level (cf. Fig. [H]). Finally, the model predicts non-zero 
anisotropic stress resulting in the difference between the perturbations of the two gravita- 
tional potentials in the conformal Newton gauge. While this effect is present already in the 
pure Einstein- aether/khronometric theory, the novel feature introduced by the aether - DM 
coupling is the persistence of the anisotropic stress with time in a wide range of momenta. 
As a consequence, its present-day power spectrum extends to wavelengths quite inside the 
current horizon (cf. Fig. [H]). 
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A qualitative comparison between the predictions of the model and those of ACDM sug- 
gests that the existing data (including local tests of gravity) have the potential to constrain 
deviations from Lorentz invariance in DM at the level of a few per cent or even better, 
Y < 0.01. This limit depends on the parameters of the aether sector and may be stronger 
or weaker for certain regions in the parameter space. Detailed numerical simulations are 
required to set the precise bounds. This work is currently in progress. 

It is worth comparing the signatures we found with the predictions of Lorentz invariant 
models where DM has unusual properties. A large class of models discussed in the literature 
[211 EHl HZl HH HHl [sol EI] includes one or several light scalar fields mediating long-range 
interactions between DM particles, referred to as fifth-forces. Depending on the model, 
the range of the fifth-force may or may not depend on the DM density. The additional 
interaction implies violation of the equivalence principle in the DM sector, similarly to our 
model. As a consequence, models with a fifth-force can also exhibit an accelerated growth of 
matter perturbations, as well as scale-dependent bias between DM and baryons. However, 
the resulting shape of the power spectrum in these models is in general different from the 
one found in this paper. The main reason is the different time dependence of the screening 
scale for the fifth-force. For example, in the simplest case of a scalar with fixed mass the 
screening momentum — the analog of our fcy (r) — is set by yU a(r) and thus grows with time. 
This means that for a given mode the fifth-force is important at early times and becomes 
screened later. On the contrary, in our case A;y(r) decreases as a(r)~^/^ and the modes enter 
into the regime of non-standard evolution at late epochs. This leads to quite different shapes 
of the transfer function and power spectrum in the two cases, compare Fig. [2] with Fig. 4 of 

One can envisage more complicated fifth-force models where the squared mass of the 
scalar field is proportional to the DM density. The screening scale in such a model will 
behave exactly as A;y(r) mimicking the effect of LV on the power spectrum. In this case 
one can try to distinguish between the two models by examining the relation between the 
anomalous growth index of perturbations and the bias factor. In the LV case this relation is 
set by Eqs. (11231) . (11261) . while in the fifth- force theories the formulas are more complicated, 
cf. Eqs. (3.11), (3.12) of [25]. The origin of this discrepancy lies in the different dynamics 
of the models. As discussed in Sec. [3l the violation of equivalence principle in our case is 
mainly due to the change in the inertial mass of the DM particles (Eq. (I30al) ). while their 
gravitational mass remains the same (Eq. (I30cl) ). The gravitational potential produced by 
DM is just the usual one. Combined with the contribution of baryons it adds up to the total 
potential that appears as the same forcing term in the hydrodynamical equations both for 
baryons and DM (Eq. fl30al) and the equivalent with F = for baryons). In other words, the 
bias between DM and baryons in the case of LV comes from the change of the inertial mass of 
DM. On the other hand, the standard fifth-force appears as an extra forcing term for the DM 
sector, without a modification of the inertial mass of DM. Thus the hydrodynamic equations 
for baryons and DM have different forcing terms with a complicated relation between them. 

Another feature of our model that distinguishes it from other fifth-force theories is the 
presence of the anisotropic stress both at early and late times. However, this effect is rather 
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small and may be outside the reach of observations (cf. Fig. [9]). 

There are several interesting questions to explore in the physics of LV DM. For large scale 
structure, it would be interesting to go beyond linear theory and understand how LV affects 
the cosmological dynamics of DM at the non-linear level. This would allow us to extend our 
results to scales shorter than 100 Mpc, where new data can be used to constrain the model. 
For this, one can resort to analytic methods in the lines of [521 |53l [511 ES] or, alternatively, use 
N-body simulations. Another open direction is related to the finding [SB] that the violation 
of equivalence principle by DM has pronounced effects on the tidal disruption of satellite 
galaxies. It seems promising to look for similar signatures at the level of galaxies and galaxy 
clusters in the model considered in this paper. One expects to encounter a rather rich 
dynamics, given that the model contains velocity-dependent interactions and automatically 
incorporates a chameleon mechanism that screens the deviations from the standard physics 
for large mass concentrations. This kind of study may lead to interesting bounds on the 
model, complementary to those discussed in the present work. 
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A Numerical procedure 

The complete system of equations for the evolution of the cosmological perturbations in 
the model consists of Einstein's equations (pTj) . equations for the khronon and DM (1931) 
and the hydrodynamical equations for the ordinary matter (1951) . Not all of these equations 
are independent. For the numerical procedure we choose two Einstein's equations (I91cp . 
(I91d|) . To these we add Eqs. (I93ap . (I93cp for the evolution of the khronon and DM velocity. 
To close the system we have to evaluate the pressure term in (I91cp . Only the radiation 
component contributes into it. The two first-order equations fl95bp governing its evolution 
can be reduced to one of second order, 

where we have introduced the density contrast. 
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Performing the Fourier decomposition, normalizing the present scale factor to one, a(ro) = 1, 
and choosing the units such that the present Hubble parameter is equal to one, Hq = 1, we 
obtain the final system of ordinary differential equations to be solved numerically: 



(3 + \ 
2 + aB 



x + 2nx + c^k\ + 



Y 



n{i-B) + n''ii + B) + 

3(2 + aB)Yndrn 



2a2 

3(2 + aB)Yndm 



X 



2a a 



V[dm] — 



2a a 

n{l + B)(l) - Btfj = 



V[dm] + ^Vydm] + + - 







^[7] + ^'^[7] + - 4^ = . 



(131a) 
(131b) 

(131c) 
(131d) 
(131e) 



A subtle point is the proper choice of the initial conditions for the system fll3ip . These 
are fixed deep inside the radiation-domination stage when the modes are superhorizon. We 
consider the initial conditions corresponding to the adiabatic mode. The latter is regular at 
r — )■ 0. Thus we write for small r: 



+ 0«r 



■,(1) 



(0) , 

17] 



[7] 



r 



Expanding Eqs. f ll3ip at r — > we obtain the relations: 



(1) 



5S = -20(r x« 



\dm\ 



(0) 



4ip 



(1) 



0(1) _ ^(1) _ 2/3x^'^ - /^(^'^"- + ^'') 



(2 + i3)0(^) + BiP^'^ - (3 + B)x 



(2) 



4./nr 



0(0) = 



(132a) 
(132b) 



(133a) 
(133b) 
(133c) 



where we have used the expansion of the scale factor at the radiation-domination epoch 
including the subleading order 

7;; , ^dm ~l~ 2 

\Ut -\ r . 



Note that Eqs. fll33ap describing the leading form of the adiabatic mode agree with those 
of Sec. 15. 1[ Additionally, the initial data must satisfy the constraint following from the (00) 
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Einstein's equation (IQlafi This gives, 



(1) (1) IWfi. (0) .(0) _^a(o) ^ 

where 5^^-^ are the constant terms in the expansion of the DM and baryon density 

contrasts. For the adiabatic mode we have: 

x(0) _ x(0) _ 3 r(0) 

Using this relation we obtain for the coefficients of the subleading terms in (11321) . 



In this way all initial conditions for the system f ll3ip are fixed in terms of the overall ampli- 
tude 0(0). 

Once the quantities entering the system (11311) are computed, we find the DM and baryon 
density contrasts by integrating the equations 

^[dm] + k'^V[dm\ - 3^" = , (135a) 
V[b] + Hvyb] -0 = 0, 5[b] + k'^V[b] -3-0 = (135b) 

with the initial conditions 

3 



By similar reasoning we also obtain the equations and initial conditions for the case of 
ACDM. As we are interested in comparing the evolution of perturbations in the proposed 
model and ACDM we can choose arbitrary normalization for 0^°), the only requirement being 
that this normalization is the same for the computations in both models. In practice we take 
0'-°'' = 1 to compute the transfer functions which we multiply by a fiat initial spectrum when 
appropriate. 

References 

[1] V. A. Kostelecky and N. Russell, Rev. Mod. Phys. 83, 11 (2011) [arXiv:0801.0287 [hep- 
ph]]. 

[2] S. Weinberg, Phys. Rev. 138, B988 (1965). 



^^It is straightforward to check that the (Oi) equation (j91bp does not produce any new constraints. 



39 



[3] S. Deser, Gen. Rel. Grav. 1, 9 (1970) [gr-qc/0411023]. 

[4] R. M. Wald, Phys. Rev. D 33, 3613 (1986). 

[5] D. Bias, J. Phys. A A 40, 6965 (2007) [hep-th/0701049]. 

[6] P. Horava, Phys. Rev. D 79, 084008 (2009) [arXiv:0901.3775 [hep-th]]. 

[7] D. Bias, O. Pujolas and S. Sibiryakov, Phys. Lett. B 688, 350 (2010) [arXiv:0912.0550 
[hep-th]]. 

[8] D. Mattingly, Living Rev. Rel. 8, 5 (2005) [gr-qc/0502097]. 

[9] T. Jacobson, S. Liberati and D. Mattingly, Annals Phys. 321, 150 (2006) [astro- 
ph/0505267]. 

[10] D. Bias, O. Pujolas and S. Sibiryakov, JHEP 0910, 029 (2009) [arXiv:0906.3046 [hep- 
th]]. 

[11] T. Jacobson and D. Mattingly Phys. Rev. D 64, 024028 (2001) [arXiv:gr-qc/0007031]. 

[12] T. Jacobson, PoS QG-PH, 020 (2007) [arXiv:0801.1547 [gr-qc]]. 

[13] B. Withers, Class. Quant. Grav. 26, 225009 (2009) [arXiv:0905.2446 [gr-qc]]. 

[14] D. Bias, O. Pujolas, S. Sibiryakov, JHEP 1104, 018 (2011). [arXiv:1007.3503 [hep-th]]. 

[15] T. Jacobson, Phys. Rev. D 81, 101502 (2010) [Erratum-ibid. D 82, 129901 (2010)] 
[arXiv: 1001. 4823 [hep-th]]. 

[16] D. Bias, O. Pujolas and S. Sibiryakov, Phys. Rev. Lett. 104 (2010) 181302 
[arXiv:0909.3525 [hep-th]]. 

[17] S. Mukohyama, Class. Quant. Grav. 27, 223101 (2010) [arXiv: 1007.5199 [hep-th]]. 

[18] D. Bias and H. Sanctuary, Phys. Rev. D 84 (2011) 064004 [arXiv:1105.5149 [gr-qc]]. 

[19] S. Groot Nibbelink and M. Pospelov, Phys. Rev. Lett. 94, 081601 (2005) [arXiv:hep- 
ph/0404271]. P. A. Bolokhov, S. G. Nibbehnk and M. Pospelov, Phys. Rev. D 72, 015013 
(2005) [arXiv:hep-ph/0505029]. 

[20] O. Pujolas and S. Sibiryakov, JHEP 1201, 062 (2012) [arXiv: 1109.4495 [hep-th]]. 

[21] M. Pospelov and Y. Shang, Phys. Rev. D 85, 105001 (2012) [arXiv: 1010.5249 [hep-th]]. 

[22] D. Bias, S. Sibiryakov, JCAP 1107, 026 (2011) [arXiv:1104.3579 [hep-th]]. 

[23] L. D. Duffy and K. van Bibber, New J. Phys. 11, 105008 (2009) [arXiv:0904.3346 [hep- 
ph]]. 

40 



[24] J. A. Frieman and B. -A. Gradwohl, Phys. Rev. Lett. 67 (1991) 2926. 

[25] B. -A. Gradwohl and J. A. Frieman, Astrophys. J. 398 (1992) 407. 

[26] J. Bovy and G. R. Farrar, Phys. Rev. Lett. 102, 101301 (2009) [arXiv:0807.3060 [hep- 
ph]]. 

[27] S. M. Carroh, S. Mantry, M. J. Ramsey- Musoh and C. W. Stubbs, Phys. Rev. Lett. 103, 
011301 (2009) [arXiv:0807.4363 [hep-ph]]; S. M. Carroll, S. Mantry and M. J. Ramsey- 
Musolf, Phys. Rev. D 81, 063507 (2010) [arXiv:0902.4461 [hep-ph]]. 

[28] J. W. Elliott, G. D. Moore and H. Stoica, JHEP 0508, 066 (2005) [hcp-ph/0505211]. 

[29] S. Dubovsky, T. Gregoire, A. Nicohs and R. Rattazzi, JHEP 0603, 025 (2006) 
[arXiv:hep-th/0512260]. 

[30] A. Adams, N. Arkani-Hamed, S. Dubovsky, A. Nicolis and R. Rattazzi, JHEP 0610, 
014 (2006) [hep-th/0602178]. 

[31] D. Bias and S. Sibiryakov, Phys. Rev. D 84 (2011) 124043 [arXiv:1110.2195 [hep-th]]. 

[32] S. M. Carroll and E. A. Lim, Phys. Rev. D 70, 123525 (2004) [arXiv:hep-th/0407149]. 

[33] B. Z. Foster and T. Jacobson, Phys. Rev. D 73 (2006) 064015 [gr-qc/0509083]. 

[34] C. M. Will, Living Rev. Rel. 9, 3 (2005) [arXiv:gr-qc/05 10072]. 

[35] C. M. Will, "Theory and experiment in gravitational physics," Cambridge, UK: Univ. 
Pr. (1993), 380 p. 

[36] B. Z. Foster, Phys. Rev. D 73, 104012 (2006) [Erratum-ibid. D 75, 129904 (2007)] 
[gr-qc/0602004]; Phys. Rev. D 76, 084033 (2007) [arXiv: 0706. 0704 [gr-qc]]. 

[37] S. Liberati and D. Mattingly, "Lorentz breaking effective field theory models for matter 
and gravity: theory and observational constraints," arXiv: 1208. 1071 [gr-qc]. 

[38] S. R. Coleman and S. L. Glashow, Phys. Rev. D 59, 116008 (1999) [hep-ph/9812418]. 

[39] I. Carruthers and T. Jacobson, Phys. Rev. D 83, 024034 (2011) [arXiv:1011.6466 [gr- 
qc]]. 

[40] J. Khoury and A. Weltman, Phys. Rev. Lett. 93, 171104 (2004) [astro-ph/0309300]. 

[41] N. Andersson and G. L. Comer, Living Rev. Rel. 10 (2005) 1 [gr-qc/0605010]. 

[42] C. Armendariz-Picon, N. F. Sierra and J. Garriga, JCAP 1007 (2010) 010 
[arXiv: 1003. 1283 [astro-ph.CO]]. 

[43] C. -P. Ma and E. Bertschinger, Astrophys. J. 455, 7 (1995) [astro-ph/9506072]. 

41 



[44] T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. 82, 451 (2010) [arXiv:0805.1726 [gr-qc]]. 

[45] G. R. Dvali, G. Gabadadze and M. Porrati, Phys. Lett. B 485, 208 (2000) [hep- 
th/0005016]. 

[46] S. Tsujikawa, Lect. Notes Phys. 800, 99 (2010) [arXiv:1101.0191 [gr-qc]]. 

[47] G. R. Farrar and P. J. E. Peebles, Astrophys. J. 604, 1 (2004) [astro-ph/0307316]. 

[48] S. S. Gubser and P. J. E. Peebles, Phys. Rev. D 70, 123510 (2004) [hep-th/0402225]; 
Phys. Rev. D 70, 123511 (2004) [hep-th/0407097]. 

[49] G. R. Farrar and R. A. Rosen, Phys. Rev. Lett. 98, 171302 (2007) [astro-ph/0610298]. 

[50] R. Bean, E. E. Flanagan, I. Laszlo and M. Trodden, Phys. Rev. D 78, 123514 (2008) 
[arXiv:0808.1105 [astro-ph]]. 

[51] J. A. Keselman, A. Nusser and P. J. E. Peebles, Phys. Rev. D 81 (2010) 063521 
[arXiv:0912.4177 [astro-ph.CO]]. 

[52] L. Amendola, Phys. Rev. D 69, 103524 (2004) [astro-ph/0311175]. 

[53] R. Scoccimarro, Phys. Rev. D 80 (2009) 104006 [arXiv:0906.4545 [astro-ph.GO]]. 

[54] F. Saracco, M. Pietroni, N. Tetradis, V. Pettorino and G. Robbers, Phys. Rev. D 82, 
023528 (2010) [arXiv:0911.5396 [astro-ph.CO]]. 

[55] P. Brax and P. Valageas, "Structure Formation in Modified Gravity Scenarios," 
arXiv:1205.6583 [astro-ph.CO]. 

[56] M. Kesden and M. Kamionkowski, Phys. Rev. D 74 (2006) 083007 [astro-ph/0608095]. 



42 



